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Abstract 

J 

Stress  near  a  crack  tip  in  plasticity  was  analyzed 
using  three  different  finite  element  modelings;  a  constant 
strain  viscoplastic  triangle,  an  eight-noded  elastic-plastic 
quadrilateral  and  a  crack  tip  elastic-plastic  singularity 
element.  The  specimen  under  consideration  was  a  center 
cracked  plate  made  from  IN-100.  The  half-crack  length  was 
0.1367  in  (3.472  mm).  An  elastic  solution  was  formulated 
and  two  different  loadings  to  generate  plasticity  were  con¬ 
sidered.  Fine  mesh  and  coarse  mesh  solutions  for  the  higher 
order  elements  were  generated  and  compared.  Comparisons 
were  made  based  on  an  equal  number  of  degrees  of  freedom  in 
two  specific  regions  referred  to  as  the  near  field  and  the 
far  field  regions. 

It  was  determined  that  the  elements  whose  elastic 
solutions  conformed  best  to  linear  elastic  fracture  mechan¬ 
ics  predictions  were  the  constant  strain  triangle  and  the 
eight-noded  quadrilateral  in  a  fine  mesh.  The  crack  tip 
element  did  not  perform  as  well  as  was  expected.  For  the 
plastic  analysis,  the  constant  strain  triangle  exhibited 
the  largest  plastic  region  and  gave  the  most  accurate 
stress  profiles.  The  eight-noded  isoparametric  element 
came  within  fifteen  percent  of  the  stress  levels  generated 
from  the  constant  strain  triangle.  The  stress  singularity 
that  is  characteristic  of  the  crack  tip  element  forced  that 
element  to  behave  unnaturally  stiff  in  plasticity. 
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Because  it  is  not  as  stiff  as  either  the  crack  tip  ele¬ 
ment  or  the  eight-noded  element,  the  constant  strain  triangle 
offered  the  most  accurate  solutions.  The  CST  solution 
required  the  least  amount  of  computational  resources. 

Though  the  isoparametric  element  mesh  was  easier  to  formu¬ 
late  and  gave  fairly  accurate  answers  elastically  and  plas¬ 
tically,  it  was  determined  that  the  constant  strain  triangle 
in  a  fine  mesh  offered  the  best  solution  to  elastic-plastic 
finite  element  problems  for  the  center  cracked  plate. 
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I .  Introduction 

The  United  States  Air  Force  has  a  vital  interest  in 
the  effects  of  cracks  on  airframe  components.  Prior  to 
1969,  the  Air  Force  determined  the  "safe  life"  of  an  air¬ 
craft  by  first  assuming  a  defect-free  structure.  Mean 
life  fatigue  analysis  was  then  used  in  which  the  projected 
lifetime  of  the  aircraft  was  divided  by  a  safety  factor. 

The  "safe  li^e"  approach  had  several  disadvantages. 
The  manufacturing  process  invariably  results  in  cracking, 
thus  invalidating  the  hypotheses  of  a  defect-free  struc¬ 
ture.  Ordinary  wear,  maintenance,  and  battle  damage  also 
create  discontinuities  in  the  structure.  The  service 
failures  of  the  F-lll  and  C-5A  emphasized  the  poor  corre¬ 
lation  between  the  safe  life  model  and  actual  aircraft 
durability.  A  new  method  of  estimating  aircraft  life  was 
needed. 

Structural  engineers  have  turned  to  Linear  Elastic 
Fracture  Mechanics  (LEFM)  to  determine  the  stresses  caused 
by  cracks  in  a  structure.  LEFM  has  worked  well;  unfortu¬ 
nately,  it  too  has  its  limitations.  LEFM  assumes  that  a 
material  will  behave  with  perfect  elasticity  until  it 
reaches  its  failure  load.  As  a  consequence,  brittle  fail¬ 
ure  will  result.  Real  materials  exhibit  plasticity  after 
the  elastic  yield  load  has  been  reached.  Higher  loads 
result  in  more  plasticity  and  greater  deviations  from  LEFM 
predictions . 
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Plasticity  occurs  around  the  tip  of  a  crack  on  a 
material  under  loading  due  to  unusually  high  stress  concen¬ 
tration  near  the  crack.  Plastic  behavior  is  nonlinear,  and 
this  results  in  irrecoverable  strains  in  the  material. 

Plastic  behavior  greatly  alters  the  material  proper¬ 
ties  in  the  region  near  the  crack  tip.  In  order  to  accu¬ 
rately  determine  the  survivability  and  durability  of  cracked 
material  under  loading,  elastic-plastic  or  elastic-visco¬ 
plastic  fracture  mechanics  must  be  utilized. 

The  primary  method  for  numerically  analyzing  stress  in 
the  vicinity  of  a  crack  is  the  "Finite  Element  Method".  To 
accomplish  this,  programs  have  been  written  incorporating 
various  finite  elements  and  the  laws  of  plastic  flow.  The 
major  purpose  of  this  thesis  is  to  determine  which  of  the 
three  different  finite  element  models  considered  is  the 
most  accurate  and  easiest  to  use. 

The  study  of  crack  growth  is  also  essential  in  the 
accurate  prediction  of  the  structural  lifetime  of  many  air¬ 
crafts.  Linear  elastic  fracture  mechanics  does  not  accu¬ 
rately  predict  crack  growth  and,  therefore,  a  plastic 
analysis  must  be  performed.  The  study  of  plastic  crack 
growth  is  beyond  the  scope  of  this  thesis.  However,  the 
information  gained  from  this  finite  element  modeling  of 
plasticity  study  will  be  useful  in  future  analyses  of 
crack  growth. 
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II .  Objective 


Plasticity  may  be  studied  in  several  ways.  A  material 
may  deform  plastically  by  the  singular  or  combined  effects 
of  load,  time  or  temperature.  A  specimen  could  be  elastic- 
perfectly  plastic,  where  the  stress  level  is  the  same  for  a 
given  plastic  strain  (see  Fig.  1),  or  it  could  exhibit 
strain  hardening,  where  stress  will  increase  with  strain 
(see  Fig.  2) .  The  model  to  be  selected  will  depend  upon 
the  desired  complexity  of  the  analysis. 

Stress  analysis  in  a  cracked  plate  is  accomplished 
today  by  the  use  of  finite  element  computer  programs. 

Since  there  are  many  finite  element  programs  available  for 
elastic-plastic  stress  evaluation,  it  would  seem  to  be  an 
easy  matter  to  select  a  program,  punch  up  data  cards,  and 
perform  the  analysis.  However,  not  all  finite  element 
programs  or  finite  element  modelings  are  alike.  The 
structural  engineer  must  determine  how  accurate  his  model 
and  program  are.  He  also  must  minimize  the  consumption  of 
resources  in  analysis  and  this  includes  computer  processor 
time,  memory  usage,  and  his  own  time  in  setting  up  the 
problem.  The  engineer  is  easily  bewildered  by  the  pro¬ 
liferation  of  programs  and  models.  As  previously  stated, 
the  primary  purpose  of  this  thesis  is  to  determine  which 
one  of  the  available  finite  element  formulations  best 
describes  a  typical  cracked  plate  problem,  in  terms  of 
accurate  results,  least  amount  of  computer  time,  and  ease 
of  use. 
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To  accomplish  this,  three  different  types  of  elements 
were  considered;  a  constant  strain  triangle  (CST) ,  an 
eight-noded  rectangular  element,  and  a  crack  tip  singular¬ 
ity  element.  Two  different  programs  were  applied  to  the 
problem;  the  VISCO  program  written  by  Captain  Terry  Hinne- 
richs,  USAF ,  and  a  program  written  by  Jalese  Ahmad  and 
V.  Papaspyropoulus  of  Battelle  Corporation,  Columbus,  Ohio, 
which  shall  be  referred  to  as  JALESE.  JALESE  utilizes  the 
8-noded  and  crack  elements  and  VISCO  incorporations  the 
constant  strain  triangle  (CST) . 

The  material  under  consideration  was  IN-100  at  1350°F 
(731. 2°C).  The  specimen  to  be  modeled  measured  1  inch 
(0.0254  m)  in  height,  1.4  in.  (0.0356  m)  in  width,  and 
0.3  in.  (0.00762  m)  in  thickness  (see  Fig.  3).  Young's 
modulus  for  IN-100  was  taken  to  be  26.3xl03  ksi  (181.2xl03 
MPa),  Poisson's  ratio  was  assumed  to  be  0.3,  and  the  yield 
stress  was  130.0  ksi  (895.9  MPa).  Because  the  specimen 
and  loading  were  symmetric,  only  the  top  quarter  of  the 
geometry  was  analyzed.  Appropriate  boundary  conditions 
were  applied  to  insure  symmetry. 

The  ratio  of  the  material  thickness,  t,  to  the  width 
of  the  specimen,  2b,  was  0.21.  This  indicated  that  the 
material  was  relatively  thin  as  compared  to  its  width. 
Therefore,  a  plane  stress  solution  was  used. 

To  guarantee  the  applicability  of  the  results,  the 
.umber  of  degrees  of  freedom  in  the  crack  tip  region  and 
the  mesh  as  a  whole  was  kept  within  a  few  percent.  For 
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comparison  purposes,  a  coarse  mesh  for  the  8-noded  and  crack 
element  in  which  the  number  of  degrees  of  freedom  in  the 
crack  tip  region  was  halved  was  also  run.  Convergence  of 
the  coarse  mesh  solution  to  the  fine  mesh  solution  was  also 
investigated . 

An  elastic  analysis  was  first  performed  on  the  speci¬ 
men,  using  the  different  finite  element  models.  A  stress 
intensity  factor  was  obtained  and  comparisons  were  made 
between  the  finite  element  modelings  and  results  from  lin¬ 
ear  elastic  fracture  mechanics.  Plastic  runs  were  then 
accomplished,  using  two  different  loadings.  Regions  of 
plasticity  and  stress  profiles  were  used  for  comparison. 
Criteria  for  program  efficiency  were  central  processor 
time,  core  memory  usage,  and  various  subjective  factors, 
such  as  ease  of  use. 


III.  Theoretical  Formulation 


Theory  of  Plasticity 

The  generalized  relationship  between  the  elastic 
stress  tensor  and  the  elastic  strain  tensor  may  be  written 
as : 


a  =  E  e  (1) 

mn  mnpr  pr 

where  o  is  the  elastic  stress  tensor,  e  is  the  elastic 
mn  pr 

strain  tensor,  and  E  are  the  elastic  constants.  For  an 

mnpr 

isotropic  material,  that  is  a  material  whose  elastic  prop¬ 
erties  are  completely  independent  of  the  orientation  of 
the  axes,  such  as  IN-100,  the  81  components  of  Emnpr  reduce 

down  to  only  two  independent  constants.  Those  constants 

3  3 

are  E,  Young's  modulus,  taken  to  be  26.3x10  ksi  (181.2x10 
MPa)  and  v,  Poisson's  ratio,  which  was  taken  to  be  0.3. 

In  compact  tensor  form,  the  equation  for  strain  versus 
stress  for  an  isotropic  material  can  be  written  as  (Ref  3) : 

e  =  ( 1+v )  o  -  v<5  a  ]  (2) 

mn  E  mn  mn  rr 

where  6  is  kronecker's  delta, 
mn 

One  would  believe  that  the  above  elastic  equation 
could  be  used  for  any  level  of  stress.  Unfortunately,  this 
is  not  true.  At  a  certain  value  of  stress,  called  the 
"yield  stress",  the  material  behavior  deviates  from  elastic 
behavior.  The  material  then  begins  to  exhibit  plasticity. 
For  IN-100,  yielding  will  occur  at  130.0  ksi  (895.9  MPa). 


For  a  uniaxial  test,  the  material  will  yield  when  the 
applied  load  equals  the  yield  stress,  and  the  specimen  will 
behave  according  to  its  uniaxial  stress-strain  curve.  The 
uniaxial  stress-strain  curve  for  IN-100  is  shown  in  Fig.  4. 
However,  for  multiaxial  stress  situations,  such  as  the 
cracked  plate  problem,  it  is  necessary  to  determine  the 
combination  of  stresses  that  will  cause  yielding.  The 
yield  criterion  indicates  when  plasticity  has  begun  (Ref  9) . 

One  of  the  most  widely  used  yield  criterion  is  the  one 
postulated  by  von  Mises  (Ref  9) .  This  criterion  states 
that  yield  in  tension  is  caused  by  maximum  distortion 
energy.  The  distortion  energy  is  given  by: 


(3) 


where  equals  distortion  energy,  G  equals  the  shear  modu¬ 
lus,  and  J2  is  the  second  invariant  of  the  deviatoric 
stress  (Ref  11) ,  defined  as 


J 


2 


=-S.  .S.  . 
2  in  in 


(4) 


where  are  the  components  of  the  deviatoric  stress 

tensor . 

This  criterion  can  also  be  viewed  as  stating  that 
yielding  will  occur  when  J2  reaches  a  critical  value,  which 
is  the  yield  stress.  For  IN-100,  the  yield  stress  in  com¬ 
pression  equals  the,  yield  stress  in  tension.  This  is  a 
characteristic  of  isotropic  materials.  Therefore,  the 
von  Mises  criterion  can  be  used  in  either  situation. 


2.0 


1.0  1.5 

:rcent  strrin 


Stress  vs .  Strain  Curve 

)0 


The  von  Mises  yield  criterion  can  be  stated  as  follows: 


F{a}  =  f(ax+ay+a2J 


+  3(t^  +t^  +t  2  )  -  (a  a  +c  a  +o  c  )  1  _  0  <5\ 

xy  yz  zx  x  y  y  z  z  x  J  o 


where  oq  is  the  uniaxial  yield  stress,  then  yielding  will 
occur  when  F(o)  >0.  The  formulation  of  the  von  Mises  yield 
criterion  is  presented  in  Mendelson  (Ref  9) . 

For  isotropic  strain  hardening,  the  yield  stress  is  a 
function  of  effective  plastic  strain  and  may  be  expressed 
as  a  segmental  linear  function  of 


) 


(6) 


where  <jq  is  the  uniaxial  yield  stress  of  the  material,  which 
is  the  stress  necessary  to  produce  yielding,  is  the  effec¬ 
tive  plastic  strain,  eQ,  and  H'  is  the  slope  of  the  uniaxial 
stress-strain  curve.  The  effective  plastic  stress  is 
defined  incrementally  by  dep  (see  Eq.  10)  . 

If  one  has  exceeded  the  uniaxial  yield  stress  of  a 
material,  the  total  strain,  emn,  can  be  expressed  as  the 
sum  of  its  plastic  and  elastic  parts. 


=  eE  +  e^3 
ran  mn 


E  p 

where  e  equals  the  elastic  strain  tensor  and  e1^  is  the 
mn  ^  mn 

plastic  strain  tensor.  A  plastic  flow  theory  which  is  con¬ 
sistent  with  the  von  Mises  criterion  is  that  the  ratio  of 


the  increment  of  plastic  strain  to  deviatoric  strain 
remains  constant,  or 


(S) 


dt? 


.A!  = 


Sij 


dA 


It  can  be  shown  (Ref  3)  that 


dA  = 


|dcp 

O'  ■ 


(9) 


where  de^  is  the  incremental  effective  plastic  strain, 
defined  as 


de^  =  yl<3£  ■  ■  de  .  . 

1  3  13  13 


(10) 


Therefore,  one  can  write: 

~  i  de 

de?.  =  |  — £  S. . 

i:  2  oe  13 

where  a  ,  called  the  effective  stress,  is  expressed  as: 
0 


(ID 


°e  = 


(12) 


Eq.  11  is  called  the  Prandtl-Reuss  equation  for  plas¬ 
tic  flow.  For  isotropic  strain  hardening,  the  Prandtl- 
Reuss  equation  becomes: 


de?  . 

13 


*3  do 
3  e 


2  H'ae  Sij 


(13) 


where  H*  is  the  slope  of  the  uniaxial  stress-strain  curve. 

This  equation  allows  one  to  solve  for  the  differential 
change  in  plastic  stress  using  known  quantities.  However, 
since  de?^  cannot  usually  be  calculated  explicitly,  an 
iterative  solution  technique  must  be  used. 
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Elastic-Plastic  Solution  Technique 


There  are  two  methods  available  for  solving  elastic- 
plastic  problems  incrementally  using  a  finite  element  com¬ 
puter  program.  The  first,  called  the  tangent  modulus 
method,  requires  a  stiffness  matrix,  [K] ,  be  constantly 
updated  to  account  for  the  effects  of  plasticity.  Since 
recalculating  K  requires  much  computer  time,  a  second 
method  for  computing  the  plastic  strain  is  desired  (Ref  13) 
The  second  method  of  calculating  plastic  strains  is 
to  divide  the  total  force  vector  into  a  plastic  load  vector 
and  an  elastic  load  vector.  This  technique,  called  the 
residual  force  method,  avoids  repeated  modifications  of 
the  stiffness  matrix.  The  residual  force  method  was  used 
in  effect  for  both  the  VISCO  and  JALESE  programs. 

When  using  the  residual  force  approach,  the  basic 
equation  of  equilibrium  for  the  finite  element  method  can 
be  expressed  as 

{F}  =  [K]  {U}  (14) 

where  { F }  is  the  generalized  force  vector,  [K]  is  the  stiff¬ 
ness  matrix,  and  { U }  is  the  generalized  displacement  vector 
and  is  broken  up  as  follows: 

tKe] {Ul1  =  {P}1  +  {Q)L  (15) 

where  K  is  the  elastic  stiffness  matrix,  P  is  the  applied 
e 

load  vector,  and  Q  is  the  effective  plastic  load  vector 
for  the  elements  in  plasticity. 


An  incremental  solution  process,  denoted  by  the  i  super¬ 


scripts,  is  used  and  there  are  two  possible  incremental 
procedures.  Both  methods  apply  an  incremental  load  and 
calculate  the  associated  values  of  stress  and  strain.  The 
initial  stress  method  approaches  equilibrium  by  iterating 
the  initial  stress  increments.  The  initial  strain  method 
approaches  compatibility  by  iterating  the  .initial  strain 
increments.  The  JALESE  program  uses  the  initial  stress 
method  to  compute  the  elastic  and  plastic  strains.  The 
VISCO  program  uses  a  time  dependent  form  of  the  initial 
strain  technique  and  is  described  in  Ref  7.  Therefore, 
the  initial  stress  method  will  now  be  examined  in  detail 
(Ref  15)  . 

The  initial  stress  method  for  the  solution  of  elas¬ 
tic-plastic  problems  is  accomplished  by  first  applying  the 
total  load  vector  and  examining  for  plasticity  by  the  yield 
criterion.  If  no  yielding  has  occurred,  one  has  a  perfectly 
elastic  problem.  If  plasticity  is  indeed  present,  then  one 
solves  for  the  elastic  and  plastic  strains  by  first  solving 
a  purely  elastic  problem  for  the  total  applied  load,  dP. 

In  this  first  step,  de1  is  determined  by  the  elastic  solu¬ 
tion.  Note  that  the  superscript  i  denotes  the  current 
increment  and  i-1  denotes  the  preceding  increment.  The 
subscript  I  denotes  the  current  iteration  and  1-1  denotes 
the  preceding  iteration. 

The  corresponding  elastic  stress,  da_T  is  calculated 
from  de .  Then,  an  initial  stress,  d0Q  is  calculated,  where 
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doo  =  daE  “  dao  (16) 

This  initial  stress  is  the  stress  required  to  maintain 

compatibility  with  the  uniaxial  stress-strain  curve.  The 

resulting  "balancing"  force  from  the  initial  stress  is 

calculated  by: 

dQ  =  J“ 3^dood  ( vol )  (17) 

dQ  is  added  to  the  applied  load,  dP,  and  the  iterative 
process  starts  over.  A  stable  solution  that  satisfies  com¬ 
patibility  and  equilibrium  is  reached  when  the  change  in 

dQ.  -  0* 

1 

Figure  5  shows  the  steps  in  the  initial  stress  algo¬ 
rithm.  Note  that  DEp  is  the  strain  to  stress  matrix  for 
elastic-plastic  solutions.  A  diagram  for  a  uniaxial  initial 
stress  solution  is  shown  in  Fig.  6. 

Viscoplastic  Material  Models 

Another  phenomenon  associated  with  problems  in  plas¬ 
ticity  is  creep.  Creep  includes  all  time-dependent  effects 
and  results  in  creep  strains,  which  are  developed  at  a 
finite  rate.  Time-independent  permanent  plastic  strains 
are  lumped  together  into  the  generalized  term  "plasticity". 
The  combination  of  time-dependent  and  time-independent 
effects  into  a  unified  plastic  flow  law  is  called  "visco¬ 
plasticity"  (Ref  14) . 

In  deriving  the  flow  law  for  viscoplasticity,  one  must 
first  decompose  the  total  strain  rate  into  its  elastic  and 
plastic  strain  rate  components. 
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For  each  load  increment  dP  in  the  plastic  range: 


1. 

ii 

tTl 

T> 

K_1(dP  +  dQ][_1) 

2. 

dex  = 

B  dUj. 

3. 

d— El  ' 

M 

(J 

'O 

Qll 

ii 

4. 

a '  = 

~  I 

2-1-1  + 

5. 

daz  = 

Sepi-i— I 

6. 

doQ  = 

4£ei  -  d£j 

7. 

SLi  = 

2i  "  d- - ! 

£-1  = 

£-1-1  +  ^1 

8. 

dQ  = 

/eT  dj<3l  d(vol) 

Process  is 

terminated  when  dQ^  =  dQ^_^ 

Figure  5.  Steps 

in  the 

Initial  Stress  Method 
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Figure  6.  Initial  Stress  Method  for  a  Uniaxial  Specimen 
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•  .  E  .  •  p 

e  .  .  =  e  .  ■  +  e1-.  • 

ID  ID  ID 


where  e. .  is  the  total  strain  rate,  e. .  is  the  elastic 
ID  ID 


strain  rate  and  e?.  is  the  plastic  strain  rate.  The  elas- 

iD 


tic  strain  rate  versus  stress  rate  equation  is  obtained 
from  the  time  derivative  of  Hooke's  law  of  elasticity.  The 
plastic  strain  rate  tensor  is  related  to  stress  rate  by  the 
time  derivative  form  of  the  Prandtl-Reuss  equation  (Ref  14) 


e?. 


ID 


=  AS  .  . 
ID 


(19) 


where  S. .  are  the  components  of  the  deviatoric  stress 


tensor.  Once  again,  incompressibility  and  isotropy  have 
been  assumed. 

It  is  not  entirely  correct  to  treat  separately  the 
plastic  strain  due  to  creep  and  the  plastic  strain  caused 
instantaneously  by  loading.  However,  the  superposition  of 
Malvern's  overstress  law  with  Norton's  law  for  secondary 
creep  has  been  proposed  as  a  unified  viscoplastic  flow  law 
(Ref  14)  into  a  finite  element  program  referred  to  as 
VISCO  (Ref  7) . 

Creep  may  be  described  as  time-dependent  deformation 
occurring  under  a  constant  strain  (Ref  13) .  The  consti¬ 
tutive  law  for  creep  is  in  the  form  of  a  strain  rate.  The 
creep  law  used  in  the  VISCO  program  is  Norton's  law  and 
may  be  expressed  as  (Ref  14) : 


de . . creep 


e^creep  = 


J -J. 


dt 


i  >B  3  ^ij 
—  Y  (o  )  — — *- 

c  e  2  o 


(20) 
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where  £?jCreep  is  the  creep  strain  rate,  e?jCreep  is  the 
creep  strain,  yc  is  the  creep  coefficient,  and  B  is  the 
creep  exponent.  Yc  and  B  are  constants  determined  from 
uniaxial  creep  test  results. 

The  Malvern  overstress  equation  relates  the  instantan¬ 
eous  plastic  strain  rate  tensor  to  stress,  and  may  be 
written  as  follows: 


"P 

£ij 


=  Yt 


_o(eP) 


-  1 


.  S.  . 
2  a 


(21] 


where  y  is  a  lluidity  constant,  a  is  the  effective  stress, 
P  ® 

and  a(ej^)  is  the  strain  hardening  yield  stress  which  has 
been  defined  in  Eq.  6.  Of  course,  this  equation  is  only 
valid  if  yield  has  been  exceeded,  that  is,  if  the  effective 
stress  is  greater  than  the  strain  hardening  yield  stress. 

If  yield  has  not  been  exceeded,  then  e?^  equals  zero. 

The  Malvern  model  can  be  adapted  for  time-independent 


elastic-plastic  solutions.  When  this  is  done,  y  may  take 
on  any  non-zero  positive  value.  The  elastic-plastic  solu¬ 
tion  is  the  steady  state  value  of  the  stresses,  strains  and 
displacements  after  the  total  load  has  been  applied.  There¬ 
fore,  if  one  multiplies  both  sides  of  Eq.  21  by  dt,  the 

3  ®ij 

resulting  expression  acting  as  the  coefficient  of  j 

e 

reduces  to  the  coefficient  for  time-independent  plasticity 
(Eq.  13) . 

VISCO  utilizes  an  Euler  linear  extrapolation  scheme  to 
integrate  the  viscoplastic  strain  rate  expressions.  The  equa¬ 
tions  for  the  Euler  extrapolation  may  be  written  as  follows: 
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(23) 


o(^)1 


°o  +  H' <£S> 1 


(24) 


The  superscript  i  refers  to  the  time  step  and  the  sub¬ 
script  i  refers  to  specific  components  of  stress  or  strain 
(Ref  7) . 

Viscoplastic  Solution  Procedure 

There  are  two  generalized  techniques  available  for 
solving  the  viscoplasticity  problem.  These  are  the  same 
two  procedures  that  were  used  to  solve  the  time-independent 
analysis.  In  the  first  method,  the  tangential  stiffness 
matrix  must  be  constantly  updated.  It  has  previously  been 
stated  that  this  method  is  inappropriate  for  use  in  a 
finite  element  computer  program.  The  second  solution  pro¬ 
cedure  uses  a  constant  elastic  stiffness  matrix  to  iter¬ 
atively  correct  the  "residual  forces"  obtained  from  the 
difference  of  elastic  and  plastic  stresses  calculated.  In 
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the  elastic-plastic  analysis,  load  was  incremented  directly 
to  satisfy  the  yield  condition  and  the  plastic  flow  law  was 
used.  In  the  elastic-viscoplastic  solution  technique,  time 
is  incremented  directly  and  load,  stress,  and  strain  are 
incremented  indirectly.  Zienkiewicz  and  Cormeau  present 
the  viscoplastic  solution  procedure,  and  it  may  be  summa¬ 
rized  as  follows: 

First,  the  time  increment,  dt1,  is  added  to  the  pre¬ 
ceding  time  t1  to  obtain  the  current  total  time,  t1. 

Then,  the  increments  of  the  plastic  strain  tensor  are  com¬ 
puted  by  the  equation: 

{d£?j}  =  {e?j}dt  (25) 

where  e. .  is  computed  from  the  Malvern  equation  (Eq.  21). 

The  plastic  strain  is  totaled  by 

{e?.}1  =  }1~1  +  { de? . } 1  (26) 

ij  ij  ij 

The  plastic  strain  rate,  c?j ,  had  been  previously  calculated 
by  the  Malvern  flow  law. 

The  plastic  load  vector,  {Q}1  ^ ,  is  computed  by 

{Q}1-1  =  f  [B]  [D]  {ef-j}1  d(vol)  (27) 

vol 

where  the  B  matrix  and  D  matrix  are  described  in  Appendix  A. 
Next,  the  current  external  load  vector  is  calculated 

by : 

{P}1  =  {P}1  dt1  +  {P}1”1  (28) 

where  P  is  the  force  r'ate  vector  (a  known  quantity)  . 
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The  nodal  displacements,  {U}1,  are  determined  from: 

{U}1  =  [K]“1({P}1  +  {Q}1"1)  (29) 

where  K  is  the  elastic  stiffness  matrix. 

The  current  total  strain  {e.^}1  is  evaluated  from  the 
strain-displacement  equation, 

{eij}1  =  [BHU)1  ’  (30) 

The  current  stress,  {a},  can  now  be  calculated  from 

{aij}1  =  [D]  ({£ij>1  ~  l^ij)1)  (31) 

The  VISCO  program  incorporates  a  variable  time  step 
procedure  to  minimize  computer  central  processor  unit  time 
usage.  This  time  step  change  is  more  fully  explained  in 
later  paragraphs.  It  is  at  this  point  in  the  program  that 
the  time  step  size  is  checked  in  terms  of  prescribed  stress 
and  strain  change  tolerances  per  time  step.  If  these  tol¬ 
erances  are  not  exceeded,  the  time  step  size  may  be 
increased  for  the  next  time  step  or  left  the  same;  if  the 
tolerances  are  exceeded,  the  time  step  size  is  reduced  and 
the  iterative  process  is  repeated  in  an  effort  to  satisfy 
the  stress  and  strain  change  tolerances.  If  the  tolerances 
have  been  satisfied,  the  iterative  process  is  repeated 
until  the  desired  simulation  time  has  been  reached. 

The  computer  algorithm  built  into  VISCO  for  estimating 
the  time  step  is  based  around  the  parameters  PQ  and  P£ 

(Ref  7),  which  are  defined  as  follows: 


(32) 
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where  atQ^  is  the  stress  tolerance,  and 
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where  eto^  is  the  strain  tolerance,  and 
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The  superscript  i  refers  to  the  time  step.  The  parameters 
Pa  and  are  calculated  for  every  element.  One  method  that 
Hinnerichs  suggested  for  changing  the  time  step  is: 


dt1  =  dt1_1/P 


(35! 


where  P  is  set  equal  to  the  largest  value  of  or  P£ . 

Another  method  that  was  suggested  for  altering  the 
time  step  that  avoids  repetitive  recalculations  uses  the 
following  equations: 


dt1  = 

0.8  dtVP 

if 

P  >  1 

dt1  = 

dt1-1 

if 

0.8  <  P  <  1 

dt1  = 

1.25  dt1-1 

if 

0.65  <  P  <  0.8 

dt1  = 

1.5  dt1”1 

if 

P  <  0.65 

Eq.  36  reduces  the  time  step  more  than  Eq.  35  if  P  is 
greater  than  1.  If  P  is  less  than  1,  Eq.  36  reduces  the 
time  step  slower  than  Eq.  35.  Hinnerichs  employed  Eq.  36 
into  the  final  form  of  VISCO. 
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Fracture  Mechanics 


Linear  Elastic  Fracture  Mechanics  (LEFM) .  Linear  elas¬ 
tic  fracture  mechanics  relates  the  stress  field  and  dis¬ 
placements  near  a  crack  tip  to  the  stress  applied  to  the 
structure  for  various  crack  geometries.  Elastic  analysis 
is  used.  The  term  used  to  describe  the  elastic  stress  and 
deformation  fields  near  a  crack  tip  is  the  stress  intensity 
factor,  K.  There  are  actually  three  different  types  of  K, 
corresponding  to  the  three  modes  of  crack  displacement; 
they  are  the  opening  mode,  sliding  mode,  and  tearing  mode 
(Ref  10)  . 

A  mode  I  (opening)  crack  occurs  when  the  stress  is 
applied  along  the  y-axis  and  the  crack  lies  along  the  x- 
axis.  The  center  cracked  plate  under  consideration  exhibits 
mode  I  behavior.  The  stress  intensity  factor,  K ^ ,  is  given 
by : 


T_  / — .  ira,  1/2 

Kt  =  a/Jia  sec  ~v-) 

I  2b 


(37) 


where  2a  is  the  crack  length  and  2b  is  the  plate  width. 

Eq.  37  is  valid  for  0.3<a/b<0.7.  Otherwise, 

_ .  i —  .2b,  ira.  1/2  ,  ..  0 , 

K  =  a/ma  —  tan  -rr-)  (38) 

I  ira  2b 

Note  that  for  the  case  of  a  =  0.1367  in.  and  b  = 

0.500  in.,  Eq.  37  is  used.  The  stresses  near  the  crack  tip, 

a  ,  a  and  x  ,  are  given  in  terms  of  KT ,  a,  radial  distance 
x'  y  xy'  y  I'  ' 

r  from  the  crack  tip,  and  angular  measure  9  from  the  crack 
tip  by  (Ref  13) : 
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The  Crack  Tip  Element .  The  crack  tip  element  accu¬ 
rately  models  the  —  singularity  that  results  from  linear 

Jr 

elastic  fracture  mechanics.  This  element  is  formed  by 

degenerating  the  eight-noded  quadrilateral  into  a  six-noded 

triangle  with  the  midside  nodes  for  the  sides  nearest  the 

crack  tip  located  at  the  quarter-chord  point  near  the  crack 

tip.  Barsoum  (Ref  2)  and  Henshell  and  Shaw  (Ref  6)  both 

describe  the  mechanism  by  which  the  —  -  singularity  is 

Jr 

formed  from  the  eight-noded  quadrilateral. 

Consider  the  quadrilateral  in  Fig.  7.  Notice  that  the 
midside  nodes  nearest  the  crack  tip  have  been  placed  at  the 
^  chord  point.  The  side  1-4  is  now  collapsed  and  the  loca¬ 
tion  of  midside  node  7  is  adjusted  such  that  it  is  at  the 
chord  position.  The  resulting  six-noded  triangle  is 
shown  in  Fig.  8.  In  this  case,  the  singularity  is  investi¬ 
gated  along  the  x-axis,  n=0. 
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Figure  7*  General  8-Noded  Quadrilateral  with 
Midside  Nodes  Moved  to  |-Chord 


J 


Figure  8.  General  Crack  Tip  Element 


The  equation  used  to  determine  x  is  given  by: 


x 


8 

E 

i=l 


N.  (5  .  ) x . 

i  in  i 


(40) 


N^  denotes  the  shape  function  for  the  eight-noded 
quadrilateral  which  is  given  in  Eq.  A-29.  By  using  the 
values  of  x  for  each  node,  the  equation  for  x  as  a  function 
of  5  for  n=0  becomes: 


x 


a 

-  i(l+5)  (l-UJ^  -  i(l+5)(l-C)i1  +  i(l-£2)-^ 

a 

+  +  |(l-C2)^  (41) 


or, 


2  1 
x  =  {^  +  2K  +  l)-~ 


(42) 


Therefore , 


5  =  -i  +  2  v! 


(43) 


Differentiating  Eq.  42  results  in 


—  =  — (i+r)  = 

35  2K  TP 


(44) 


where  the  value  of  5  in  terms  of  x  given  by  Eq.  43  has  been 

3x 

substituted.  Notice  that  x-=-  equals  zero  at  the  crack  tip 

° 


(x=0) .  Now  . 


=  l1 


(45) 


where 


detJ 


3 (x,y) 

3 (5,n) 


(46) 
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Therefore,  Eq.  44  makes  the  Jacobian  singular  at  x=0, 
5=-l.  Considering  only  the  displacements  at  nodes  1  and  6, 


the  displacement  u  along  the  line  1-6  is  given  by: 

u  =  J(C2-l)u1  +  |(l+£2)u6 


(47) 


The  strain  in  the  x  direction  then  becomes: 


3u  _  T-l  3u  =  3£  au 
£x  3x  =  3 £  3x  35 


=  (-1  Vi +  l)“i  *  (-  vf +  2)u 


6 

(48) 


Therefore,  a  —  strain  singularity  along  the  line  1-2 
/r 

has  been  created  by  the  use  of  the  crack-tip  singularity 

element.  Since  stress  is  proportional  to  strain,  a  — 

/r 

stress  singularity  also  has  been  created.  Since  the 

Jacobian  always  goes  to  zero  at  the  crack  tip,  one  can 

choose  any  value  of  theta  and  there  will  be  a  —  singular- 

/r 


ity  for  stress  and  strain  along  that  r. 

If  the  crack  tip  is  surrounded  by  singularity  elements 

(see  Fig.  9) ,  then  one  can  choose  any  value  of  theta  and 

have  a  —  singularity  for  stress.  Therefore,  surrounding 
/r 

the  crack  tip  with  singularity  elements  effectively  simu¬ 
lates  the  —  singularity  from  the  theory  of  linear  elastic 
/r 

fracture  mechanics. 

Plasticity  within  the  Singularity  Elements .  Any  finite 
element  will  depict  plastic  action  when  the  effective  stress 
within  the  element  has  reached  the  yield  stress.  When  this 
occurs,  the  stresses  and  strains  within  the  element  can  be 
determined  only  through  the  iterative  plastic  analysis 
described  previously.  Linear  elastic  fracture  mechanics 


X 

NCrack  Tip 


Figure  9«  Crack  Tip  Surrounded  with  Crack 
Tip  Elements 


can  no  longer  be  used  to  determine  the  stres.:"^  near  the 

crack  tip.  In  other  words,  the  total  stress  will  not 

increase  as  a  function  of  —  as  one  approaches  the  crack 

/r 

tip.  Eventually,  the  effective  stress  for  the  elements 

near  the  crack  tip  will  reach  the  maximum  value  that  is 

permitted  for  the  material  under  plasticity,  which  for 

IN-100  is  taken  as  164.0  ksi  (1130.1  MPa).. 

The  crack  tip  singularity  element  performs  irregularly 

in  the  crack  tip  region  when  the  yield  stress  is  reached  due 

to  the  fact  that  the  element  is  being  forced  to  obey  the 

—  singularity  rule  even  when  it  no  longer  applies.  There- 
/r 

fore,  spurious  values  of  the  stresses  and  strains  for  the 
crack  tip  element  in  plasticity  should  be  seen  in  the  finite 
element  computer  analysis. 
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IV.  Results  and  Discussion 

Time  Independence  vs .  Time  Dependence 

The  JALESE  finite  element  program  performs  a  time- 
independent,  elastic-plastic  analysis  and  uses  the  Prandtl- 
Reuss  equations  to  determine  plastic  strain.  The  VISCO 
finite  element  program,  on  the  other  hand,  performs  a  time- 
dependent,  elastic-viscoplastic  analysis  and  uses  a  time- 
dependent  equation  to  determine  plastic  strain  such  as  the 
Malvern  flow  law  and  Norton's  law  for  creep.  When  solving 
plasticity  problems  using  JALESE,  one  first  had  to  input  a 
small  elastic  load  and  the  program  then  multiplied  the  load 
until  one  element  became  plastic.  Then  it  was  necessary  to 
input  small  increments  of  load  in  order  to  arrive  at  a 
plastic  solution.  When  using  the  VISCO  program,  one  had 
to  input  the  total  desired  load  and  a  rate  at  which  the 
load  was  to  be  applied.  The  program  used  a  time-stepping 
procedure  to  increment  the  load  and  compute  the  plastic 
strain  (see  the  section  on  Viscoplastic  Solution  Procedure 
in  the  Theoretical  Formulation) . 

The  time-dependent  solution  procedure  used  in  VISCO 
was  adapted  for  time-independent  problems  by  first  setting 
the  creep  coefficient  equal  to  zero  (Ref  8)  .  Next,  7^  was 
assigned  a  non-zero  positive  value  which  produced  a  steady 
state  solution.  Finally,  the  load  was  applied  quickly; 
therefore,  a  high  force  rate,  such  as  20%  of  the  load  per 
second,  was  used. 
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Mesh  Arrangement 

The  first  objective  of  the  present  research  was  to 
verify  that  the  VISCO  and  JALESE  programs  were  working. 

This  was  accomplished  .by'  testing  a  uniaxial  specimen  using 
a  minimum  number  of  elements.  If  the  resulting  plot  of 
stress  versus  strain  matched  the  uniaxial  curve,  the  pro¬ 
gram  was  determined  to  be  verified.  For  the  constant 
strain  triangle  in  the  VISCO  program,  only  two  elements 
were  used  (see  Fig.  10) .  The  JALESE  program  was  tested 
with  two  three-noded  triangles  and  one  eight-noded  quadri¬ 
lateral  (see  Fig.  11) .  The  triangle  available  in  the 
JALESE  program  is  not  a  constant  strain  triangle  but  a 
four-noded  quadrilateral  condensed  down  to  a  triangle. 

It  uses  the  shape  functions  and  Gaussian  integration 
schemes  that  a  quadrilateral  uses  and  not  the  simple  func¬ 
tions  and  integration  that  the  CST  of  the  VISCO  program 
uses . 

For  each  of  the  test  cases,  the  effective  stress  was 
plotted  against  the  effective  strain.  The  results  were 
found  to  correspond  very  closely  to  the  uniaxial  stress- 
strain  curve  for  IN-100.  Since  this  satisfied  the  verifi¬ 
cation  criterion,  the  programs  were  held  to  be  in  working 
condition . 

The  next  task  that  had  to  be  accomplished  before  the 
comparison  tests  could  be  run  was  to  devise  a  suitable 
finite  element  mesh.  For  the  constant  strain  triangle,  the 
mesh  that  was  used  was  identical  to  the  one  presented  by 
Hinnerichs  (Ref  7) .  This  mesh  will  be  discussed  in  later 
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paragraphs.  For  the  8-noded  and  crack  tip  elements,  it  was 
necessary  to  slowly  develop  a  finite  element  mesh  by  first 
devising  a  very  coarse  mesh  using  only  95  nodes,  170  degrees 
of  freedom  and  26  elements  for  the  crack  tip  element  (see 
Fig.  12) ,  and  107  nodes,  192  degrees  of  freedom  and  28  ele¬ 
ments  for  the  8-noded  element  (see  Fig.  13) .  The  elastic 
stress  intensity  factor  for  each  of  these  cases  was  deter¬ 
mined,  and  the  results  were  compared  to  the  stress  intensity 
factor  for  the  constant  strain  triangle,  fine  mesh  given  in 
Hinnerichs  as  shown  subsequently  (Ref  7) .  The  results  were 
also  compared  with  the  theoretical  stress  intensity  factor 
obtained  from  linear  elastic  fracture  mechanics. 

Once  it  had  been  verified  that  the  very  coarse  mesh  pat¬ 
terns  were  giving  fairly  accurate  linear  elastic  answers,  it 
then  became  necessary  to  devise  the  refined  mesh  patterns  for 
the  crack  tip  element  mesh  and  the  8-noded  element  mesh.  Bar- 
soum  (Ref  2)  states  that  a  mesh  pattern  that  propagates  radially 
from  the  crack  tip  works  the  best  for  elastic  analysis  because 
a  radial  pattern  allows  one  to  easily  change  element  size. 

Such  a  radial  pattern  was  used  in  the  fine  meshes.  The  same 
mesh  pattern  was  used  for  elastic  and  plastic  analysis  since 
the  mesh  could  not  be  altered  in  the  middle  of  the  program 
run.  The  mesh  patterns  are  illustrated  in  Figs.  14-16.  The 
crack  tip  region  (§•  <  0.005)  meshes  are  shown  in  Figs.  17-19. 

Cl 

Since  the  results  of  the  eight-noded  and  crack  tip  ele¬ 
ment  runs  were  to  be  compared  with  results  from  the  constant 
strain  triangle,  it  was  necessary  to  keep  the  number 
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Figure  12.  Crack  Tip  Element  Test  Mesh 
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Figure  16.  CST  Fine  Mesh 
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of  degrees  of  freedom  in  the  various  regions  the  same.  By 
doing  this,  one  is  able  to  compare  each  type  of  element  to 
their  primary  characteristics.  The  number  of  degrees  of 
freedom  is  the  important  parameter  because  the  number  of 
displacement  equations  to  be  solved  equals  the  number  of 
degrees  of  freedom.  Equal  degrees  of  freedom  in  given 
regions  isolates  the  effect  of  element  types  as  discussed 
subsequently.  Therefore,  for  comparable  results,  it  is 
important  that  the  number  of  degrees  of  freedom  in  each 
test  case  be  kept  the  same. 

Table  I  shows  the  number  of  degrees  of  freedom  in 
each  region.  Notice  that  overall,  the  crack  tip  element 
mesh  had  60  elements,  205  nodes  and  382  degrees  of  freedom, 
the  eight-noded  element  mesh  had  58  elements,  207  nodes, 
and  382  degrees  of  freedom,  and  the  constant  strain  tri¬ 
angle  mesh  had  355  elements,  211  nodes,  and  378  degrees  of 
freedom.  The  crack  tip  element  mesh  and  the  eight-noded 
element  mesh  had  exactly  the  same  number  of  degrees  of 
freedom  and  the  constant  strain  triangle  mesh  had  only  1% 
fewer  degrees  of  freedom. 

The  two  most  important  regions  for  which  the  number  of 
degrees  of  freedom  must  be  under  close  scrutiny  are  the  far 
field  region,  for  this  is  where  the  stresses  are  applied, 
and  the  crack  tip  region,  since  this  is  where  the  most 
important  stress  analysis  will  occur.  In  the  crack  tip 
region,  the  crack  tip  element  mesh  and  the  eight-noded 
element  mesh  had  exactly  the  same  number  of  degrees  of 
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TABLE  I 


Regions  and  Degrees  of  Freedom,  Fine  Meshes 


region 

see  sketch  below 

r/a  (approx . ) 

no.  of  d.c 

).  f  .  (approx.  ) 

crack  el. 

8-noded 

CST 

crack  tip  (CT) 

0.05 

153 

153 

151 

near  field  (NF) 

0. 05<r/a<0. 5 

113 

113 

125 

mid  field  (MF) 

0 . 5<r/a<l . 0 

57 

57 

37 

transition  field 
(TF) 

1 . 0<r/a<2 . 0 

8 

8 

12 

far  field  (FF) 

2 . 0<r/a<4 . 0 

53 

53 

59 

Totals 

_ _ _ 

fine  mesh  type 

no.  of 
elements 

no.  of 
nodes 

mmm 

1 

crack  tip  element 

60 

205 

382 

8-noded  element 

58 

207 

382 

CST 

355 

211 

378 

(jT 


freedom,  and  the  constant  strain  triangle  mesh  had  only  4% 
more  degrees  of  freedom.  In  the  far  field  region,  the  crack 
tip  element  mesh  and  the  eight-noded  element  mesh  had  the 
same  number  of  degrees  of  freedom,  and  the  constant  strain 
triangle  mesh  had  only  10%  more  degrees  of  freedom.  Large 
differences  occur  in  the  intermediate  regions.  This  is 
unavoidable.  The  reductions  in  element  size  for  the  eight- 
noded  and  crack  tip  element  meshes  occurred  differently 
than  for  the  constant  strain  triangle  mesh. 

After  the  tests  with  the  fine  meshes  were  completed,  a 
coarse  mesh  in  which  the  number  of  degrees  of  freedom  in 
the  crack  tip  region  was  halved  was  tried  for  the  crack  tip 
element  and  the  eight-noded  element  meshes.  The  crack  tip 
region  arrangement  for  the  coarse  meshes  are  shown  in  Fig. 

20  and  21.  The  coarse  meshes  were  used  to  validate  the 
convergence  of  fine  mesh  solutions. 

Elastic  Analysis-Stress  Intensity  Factor 

The  stress  intensity  factor  describes  the  magnitude  of 
the  elastic  stress  field  in  the  crack  tip  region  (Ref  10) . 

was  defined  in  Eq.  37.  The  relation  between  and 

stress  for  any  r  or  angle  0  is  given  in  Eq.  39. 

Chan  (Ref  4)  gives  a  nondimensionalized  form  of  the 
elastic  stress  intensity  factor.  This  may  be  written  as: 

KI 

Kt  =  -  (49) 

For  the  center-cracked  plate  with  a=0.1367  in.  and 
b=0.5  in.,  Kj  is  equal  to  1.86. 


Kj  was  calculated  numer- 


Numerical  Evaluation  of  K ^ . 
ically  from  the  elastic  finite  element  solutions.  A  radial 
line  was  drawn  from  the  crack  tip  at  a  small  angle,  which 
was  taken  to  be  7.5°.  In  order  to  determine  ,  an  elastic 
analysis  was  first  performed  with  the  finite  element  prob¬ 
lem.  Values  of  a  were  taken  at  several  stations  of  r  in 

y 

the  crack  tip  region.  For  each  value  of  a  ,  was  deter¬ 
mined  by  the  following  equation,  which  is  a  form  of  Eq.  39: 

Kj  =  0^  /2irr  [cos  j(l  +  sin  ~  sin  ^-)  ]  ^  (50) 

The  value  of  obtained  was  nondimensionalized  by 
dividing  it  by  (the  applied  stress) .  The  resulting 
nondimensionalized  stress  intensity  factor,  KI#  was  then 
plotted  against  its  value  of  r/a.  This  was  done  for  all  of 
the  points  in  the  crack  tip  region.  The  resulting  curve  of 
vs.  r/a  was  extrapolated  to  the  K^.  axis  (r/a=0)  .  This 
was  the  elastic  value  of  the  nondimensionalized  stress 
intensity  factor  for  the  problem. 

The  graphs  of  nondimensionalized  stress  intensity  fac¬ 
tor  versus  r/a  are  presented  in  Fig.  22  and  Fig.  23.  Table 
II  presents  the  extrapolated  value  of  the  nondimensionalized 
stress  intensity  factor  for  each  mesh. 

From  Table  II,  it  is  shown  that  the  constant  strain 
triangle  and  the  eight-noded  fine  mesh  both  generate  the 
same  value  of  Kj  that  was  calculated  using  linear  elastic 
fracture  mechanics.  The  eight-noded  element  in  a  fine  mesh 
and  the  constant  strain  triangle  in  a  fine  mesh,  therefore, 
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Figure.  22.  Elastic  Analysis  (Kj),  Coarse  Meshes 
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C Crack  Element 


□8-Noded  Element 
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RADIAL  DISTANCE  (R/A) 


Figure  23.  Elastic  Analysis  (KT),  Fine  Meshes 


TABLE  II 

Elastic  Analysis,  Values  of  K ^  for  Fine  and  Coarse  Meshes 


type  of  element  and  mesh 

KI 

crack  tip  element  -  coarse 

2.34 

8-noded  element  -  coarse 

1.95 

crack  tip  element  -  fine 

1.95 

8-noded  element  -  fine 

1.86 

CST  -  fine 

1.86 

T 


were  the  best  elements  for  solving  the  elastic  problem.  It 
is  noted  from  Fig.  22  and  Fig.  23  that  the  crack  tip  ele¬ 
ment  does  effectively  simulate  the  —  singularity,  K 

/? 

becomes  very  large  near  the  crack  tip.  However,  the 
extrapolated  values  of  for  this  element  were  not  as 
close  to  the  K  values  for  the  constant  strain  triangle  and 
the  eight-noded  element.  The  crack  tip  singularity  element 
did  not  perform  as  well  elastically  as  Barsoum  predicted 
(Ref  2) .  This  was  true  both  for  the  fine  and  coarse 
meshes . 

There  was  a  20%  difference  between  the  values  of  K^. 
for  the  crack  element  fine  mesh  and  the  crack  element 
coarse  mesh.  The  K  values  for  the  eight-noded  fine  and 
coarse  meshes  also  differed  by  20%.  Therefore,  the  coarse 
mesh  elastic  solutions  converged  to  the  fine  mesh  elastic 
solutions . 

Plastic  Analysis 

Two  different  values  of  loading  were  used  for  the 
plastic  analysis;  10.896  Klb  (48.465  kN)  and  16.060  Klb 
(71.435  kN) .  The  applied  stresses  therefore  were  36.320  ksi 
(250.29  MPa)  and  53.533  ksi  (368.91  MPa).  These  loadings 
corresponded  .to  values  of  K  of  25.0  ksi /In  (27.4  MPa/m)  and 
36.8  ksi/In  (40.41  MPa^in)  .  As  stated  in  Appendix  B,  the 
loadings  had  to  be  applied  as  equivalent  nodal  forces. 

This  procedure  is  outlined  in  Appendix  B. 

The  comparison  criteria  for  the  plastic  analysis  were 
plastic  regions,  profiles  of  vs.  r/a  for  several  angles. 
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Figure  24.  Stress  Profiles,  0=  7-5°> 
P  =  10.896  klb 
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Figure  25.  Stress  Profiles,  0=  7-5°. 
P  =  16.060  klb 
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Figure  27.  Stress  Profiles,  0=  45° , 
P  =  16.060  klb 
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Figure  28.  Stress  Profiles,  0=  75°, 
P  =  IO.896  klb 


5: 


0  Crack  Element 


"  i  r  i - 1 - 1 - 1 - 1 - 1 - 1 - i - 1 - j - r  i - 1 - 1 - 

.0  0.4  0.8  1.2  1.6 

DISTANCE  FROM  CRACK  TIP  (R/A) 


Figure  29.  Stress  Profiles,  9=  ?5°  , 
P  =  16.060  klb 
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Figure  30.  Stress  Profiles,  @=  112.5°, 
P  =  10.896  kit 


EFFECTIVE  STRESS  (KSI) 


DISTANCE  FROM  CRACK  TIP  (R/A) 


Figure  .  Stress  Profiles)  f?=112.5°> 
P  =  16.060  klb 
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Figure  33-  Stress  Profiles,  0=  157-5  » 
P  =  16.060  klb 


and  effective  stress  profiles.  An  element  was  determined 
to  have  gone  plastic  if  its  effective  stresses  had  exceeded 
the  yield  stress.  This  is  simply  a  restatement  of  the  von 
Mises  yield  criterion. 

The  stress  profiles  shown  as  oq  vs.  r/a  at  various 
angles  9  are  shown  in  Fig.  24  through  Fig.  33.  The  plastic 
regions  are  shown  in  Fig.  34  through  Fig.  39. 

The  diagrams  of  the  plastic  region  show  that  the  con¬ 
stant  strain  triangle  mesh  had  the  largest  region  of  plas¬ 
ticity  at  both  load  levels.  The  plastic  regions  for  the 
other  elements  covered  roughly  the  same  areas.  It  is  not 
surprising  that  the  CST  mesh  would  show  the  greatest  amount 
of  plasticity;  the  constant  strain  triangle  is  more  flexible 
than  either  the  crack-tip  singularity  element  or  the  eight- 
noded  quadrilateral. 

The  stress  profiles  show  that  the  crack-tip  singular¬ 
ity  element  generated  extraordinarily  high  stresses  at  the 
crack  tip.  Figure  40  and  Fig.  41  show  the  curves  of  x- 
displacement  vs.  distance  from  the  crack  tip  for  each  type 
of  element  at  the  two  different  loadings.  The  slope  of 
these  curves  near  the  crack  tip  was  determined,  and  this 
value  was  divided  by  the  crack  length  to  obtain  3y/3x, 
which  is  c  ,  or  the  strain  in  the  x-direction.  The  results 
were  compiled  in  Table  III.  It  is  noted  that  the  crack-tip 
singularity  element  exhibited  strains  that  were  an  order  of 
magnitude  greater  than  either  the  eight-noded  quadrilateral 
or  the  constant-strain  triangle.  Since  the  singularity 


TABLE  III 


Strains  Near  the  Crack  Tip,  0=0°,  Fine  Meshes 


load  (klb) 

fine  mesh  type 

near  crack  tip 

10.896 

crack  tip  element 

40 . 0xl0~3 

-3 

10.896 

8-noded  element 

9.52x10 

10.896 

CST 

3 . 55xl0~3 

16.060 

crack  tip  element 

60 . 92xl0-3 

16.060 

8-noded  element 

7 . 7  5xl0_3 

16.060 

CST 

3 . 98xl0~3 
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0.04  0.08  0.12  0.18 


DISTANCE  FROM  CRACK  TIP  (R/A) 


Figure  4l  .  Displacement  Profiles,  9=  0° , 
p  -  16.060  klb 


that  is  characteristic  of  the  crack  tip  element  creates  high 

levels  of  stress  near  the  crack  tip,  and  since  stress  is 

related  to  strain,  the  —  singularity  forces  the  crack  tip 

/r 

element  to  behave  with  unnatural  stiffness  in  the  plastic 
regime. 

It  also  should  be  noted  that  the  eight-noded  element 
produced  strains  that  were  about  three  times  greater  than 
the  strains  in  the  CST  near  the  crack  tip.  This  is  due  to 
the  fact  that  higher  order  elements  are  naturally  stiffer 
than  the  constant  strain  triangle. 

The  stress  profiles  show  that  though  the  eight-noded 
element  produced  stresses  that  were  higher  than  the  con¬ 
stant  strain  triangle,  it  was  not  as  far  off  as  was  the 
crack  tip  element.  The  eight-noded  element  gave  stresses 
that  were  within  15%  of  the  stress  values  of  the  CST, 
except  in  the  elements  nearest  the  crack  tip.  This  further 
demonstrates  that  the  eight-noded  element  was  stiffer  than 
the  constant  strain  triangle,  but  it  was  not  as  stiff  as 
the  singularity  element. 

Each  of  the  finite  element  modelings  accurately  dis¬ 
played  the  stresses  near  the  stress-free  boundary  at  0=180°. 
At  0=157.5°  (Fig.  32  and  Fig.  33),  very  high  stresses  were 
shown  in  the  element  nearest  the  crack  tip.  The  stresses 
rapidly  dropped  off  away  from  the  crack  tip.  In  the  region 
of  0.1  <  r/a  <  0.6,  the  eight-noded  element  meshes  showed 
markedly  higher  stresses  than  the  CST  or  crack  element 
meshes.  This  occurred  because  there  were  two  elements  near 
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CFine  Mesh 


1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - r 

0.8  1.2 
RADIAL  DISTANCE  (R/A) 


I — 1 — 1 
1.6 


Figure  42 .  Convergence  of  Coarse  to  Fine  Mesh 
Crack  Tip  Element  Mesh 


EFFECTIVE  STRESS  (KSI) 


the  crack  tip  for  the  crack  element  mesh  and  several  ele¬ 
ments  near  the  crack  tip  for  the  CST  mesh,  while  there  was 
only  one  element  near  the  crack  tip  for  the  eight-noded 
mesh.  Therefore,  it  was  easier  for  the  CST  and  crack  ele¬ 
ment  meshes  to  adjust  to  the  stress-free  boundary  than  the 
eight-noded  element  mash. 

Figure  42  and  Fig.  43  show  the  convergence  of  the 
coarse  meshes  to  the  fine  meshes.  For  r/a  >  0.70,  the 
coarse  mesh  and  the  fine  mesh  for  both  types  of  elements 
each  gave  identical  results.  However,  for  r/a  <  0.7,  the 
coarse  meshes  gave  stresses  that  ranged  from  three  to  seven 
percent  higher  than  the  fine  meshes. 

The  diagrams  of  the  plastic  region  show  that  the  con¬ 
stant  strain  triangle  produced  the  largest  area  of  plastic¬ 
ity.  Since  the  number  of  degrees  of  freedom  in  each  mesh 
was  kept  the  same  and  the  number  of  degrees  of  freedom  in 
the  crack  tip  region  was  kept  the  same,  the  larger  region 
of  plasticity  in  the  constant  strain  triangle  mesh  was 
caused  by  the  fact  that  the  constant  strain  triangle  was  not 
as  stiff  as  either  the  eight-noded  element  or  the  crack  tip 
element,  and  not  because  one  mesh  was  inherently  more 
restrained  than  another.  Because  of  its  close  agreement 
elastically  with  linear  elastic  fracture  mechanics  and 
because  of  its  lower  stiffness,  the  constant  strain  tri¬ 
angle  is  preferred  over  the  eight-noded  element  and  the 
crack  tip  element  for  modeling  elastic-plastic  problems. 
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The  higher-order  elements  do  have  some  advantages,  how¬ 
ever.  Since  a  mesh  involving  the  higher-order  elements  has 
fewer  elements  than  the  constant  strain  triangle  mesh,  it 
requires  less  work  on  the  part  of  the  user  to  input  the  mesh 
for  the  higher-order  elements.  The  eight-noded  rectangular 
program  also  requires  that  only  corner  nodes  be  inputed  while 
the  constant  strain  program  requires  that  the  user  supply  the 
coordinates  of  all  of  the  nodes.  Since  less  cards  are  punched 
for  the  higher-order  element  mesh  than  for  the  constant  strain 
triangle  mesh,  there  is  less  chance  of  error;  and  since  there 
are  fewer  cards  for  the  isoparametric  elements,  it  is  easier 
to  check  the  mesh.  Overall,  one  would  produce  results  more 
quickly  by  using  higher-order  elements  than  by  using  the  con¬ 
stant  strain  triangle  in  terms  of  elapsed  time  from  the 
start  to  the  finish  of  the  problem.  However,  it  appears  the 
elastic-plastic  resulting  stresses  will  be  more  realistic  if 
one  uses  the  constant  strain  triangular  element  rather  than 
the  eight-noded  element  or  crack-tip  singularity  element. 

Computational  Comparisons 

The  last  comparison  to  be  made  involves  the  difference 
in  computer  resources  required  between  the  two  programs.  A 
test  case  was.  set  up.  A  uniaxial  specimen  was  divided 
into  two  triangular  elements  acting  under  a  uniaxial  load 
of  134.0  ksi  (923.4  MPa).  The  JALESE  and  VISCO  programs 
were  then  used  for  comparison.  Triangular  elements  were 
chosen  because  both  programs  incorporated  triangular  ele¬ 
ments.  However,  the  triangular  elements  used  in  the  JALESE 
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program  were  not  constant  strain  triangles  but  degenerated 
four-noded  rectangles.  Once  again,  the  creep  coefficient 
was  set  equal  to  zero  and  the  load  was  applied  very  quickly. 
The  results  are  shown  in  Table  IV. 

The  JALESE  program  used  30%  more  central  processor 
time  than  the  VISCO  program  in  the  test  cases.  This 
occurred  because  the  JALESE  program  utilizes  Gaussian  quad¬ 
rature  to  perform  the  integration  necessary  to  obtain  the 
stiffness  matrix  while  the  VISCO  program,  which  incorporates 
a  constant  strain  triangle,  can  evaluate  the  integrals 
exactly  since  only  constants  are  integrated.  The  JALESE 
program  also  required  32300  (octal)  more  core  memory  to 
load  than  the  VISCO  program.  This  is  because  the  VISCO 
program  stores  the  K  matrix  in  a  more  efficient  manner. 

The  greater  usage  of  core  memory  and  central  processor  time 
slowed  down  the  turnaround  time  of  JALESE  as  compared  to 
VISCO.  It  also  makes  the  JALESE  runs  more  expensive  than 
the  VISCO  runs. 

If  the  default  value  of  stress  and  strain  tolerance 
(0.001)  was  used  with  the  crack  tip  element,  fine  mesh,  in 
the  JALESE  program,  the  iterative  initial  stress  solution 
for  the  plastic  strains  would  not  converge.  The  stress 
tolerance  had  to  be  increased  to  0.005  in  order  to  obtain 
a  solution.  Since  the  program  only  printed  out  four  sig¬ 
nificant  figures  for  the  stress  and  strain  solution,  this 
change  of  the  tolerance  did  not  make  a  significant  differ¬ 
ence  . 
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TABLE  IV 

Computational  Resource  Usage  Comparisons, 
Two  Triangular  Elements,  Uniaxial  Test, 
JALESE  and  VISCO 


Program  Name 

CM  required 

CP  usage 

10  time 

Cost 

(octal  words) 

(seconds) 

(seconds) 

(dollars) 

JALESE 

205100 

5.560 

25.424 

1.04 

VISCO 

152600 

4.273 

25.678 

0.94 
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It  also  should  be  noted  that  though  it  was  easier  to 
input  data  into  the  JALESE  program,  it  was  harder  to  obtain 
the  desired  levels  of  stress  using  JALESE  over  VISCO. 

JALESE  required  that  the  user  supply  incremental  units  of 
loading  after  plasticity  has  been  achieved.  When  using 
VISCO,  the  total  desired  load  was  inputed  in  one  data  card. 
VISCO  used  time  increments  to  increment  the  load. 
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V.  Conclusions 


It  appears  that  for  at  least  the  center  cracked  plate 
the  best  element  to  use  for  elastic-plastic  analysis  was 
the  constant  strain  triangle.  Elastically,  the  constant 
strain  triangle  gave  results  that  agreed  exactly  with  lin¬ 
ear  elastic  fracture  mechanics.  The  eight-noded  quadri¬ 
lateral  also  performed  well  elastically;  its  results  also 
agreed  with  linear  elastic  fracture  mechanics  predictions. 
The  crack  tip  singularity  element,  however,  did  not  give 
accurate  elastic  results  despite  its  unique  ability  to  sim 
ulate  the  singularity  that  occurs  in  the  crack  tip  region. 

Plastically,  the  constant  strain  triangle  performed 
much  better  than  either  the  eight-noded  quadrilateral  or 
the  crack  tip  singularity  element.  The  constant  strain 
triangle,  in  a  mesh  with  the  same  number  of  degrees  of 
freedom  as  the  other  elements,  exhibited  larger  plastic 
regions  and  more  realistic  stresses  near  the  crack  tip. 
This  was  verified  by  Hinnerichs  by  comparisons  made  with 
experimental  measurements  (Ref  7) .  The  eight-noded  quadri 
lateral  gave  erroneously  high  stresses  near  the  crack  tip 
because  of  its  inherent  stiffness.  The  crack-tip  singu¬ 
larity  element  gave  even  higher  stresses  because  of  its 
imposed  singularity.  This  singularity  made  the  crack  tip 
elements  unnaturally  stiff.  They  did  not  model  plasticity 
very  well. 

Though  the  eight-noded  quadrilateral  element  and  the 
crack-tip  singularity  element  meshes  were  easier- to  input 
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and  check  than  the  constant  strain  triangle  element  mesh, 
it  was  easier  to  obtain  the  desired  stress  level  using  the 
triangle  because  the  program  that  was  used  for  the  triangu¬ 
lar  element  incremented  with  time  viscoplastically  and  not 
with  user-supplied  load  increments.  The  VISCO  program  that 
incorporated  the  constant  strain  triangle  viscoplastically 
required  less  core  memory  space  and  less  computer  process 
or  time  than  the  JALESE  program  that  contained  the  higher- 
order  element  approach. 

The  constant  strain  triangle  also  modeled  the  stress- 
free  region  better  than  either  the  eight-noded  element  or 
the  crack  tip  element.  The  constant  strain  triangle  with 
its  great  flexibility  material-wise  provided  a  very  accu¬ 
rate  model  of  the  elastic  and  plastic  strains  near  the 
crack  tip. 

One  recommendation  for  future  work  is  that  a  time 
dependency  be  incorporated  into  the  JALESE  program  and  that 
the  comparisons  be  made  viscoplastically.  Creep  should  be 
studied  in  addition  to  the  instantaneous  plastic  strains. 

A  constant  strain  triangle  should  also  be  incorporated  into 
the  JALESE  program  so  that  the  comparisons  could  all  be 
made  with  one  computer  program.  Another  suggested  future 
project  is  that  crack  growth  comparisons  should  be  made. 

In  summary,  the  eight-noded  and  crack-tip  singularity 
elements  presented  in  JALESE  perform  well  elastically,  with 
the  crack  element  being  slightly  less  accurate.  The  plas¬ 
tic  results  of  the  crack  element  and  the  eight-noded 
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element  are  disappointing.  The  constant  strain  triangle, 
with  its  great  simplicity  for  understanding,  provides 
extremely  accurate  answers  in  the  elastic  and  plastic 
regimes.  Though  some  extra  work  is  required  to  set  up  the 
mesh,  it  is  well  worth  it  in  order  to  obtain  comparable 
accuracy  with  faster  computer  turnaround  time. 
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Appendix  A:  The  Finite  Element  Method 


Introduction 

The  basic  philosophy  of  the  finite  element  is  that  an 
approximate  solution  for  a  problem  too  complex  to  be  solved 
exactly  can  be  found  by  dividing  the  field  of  interest  into 
a  number  of  discrete  regions  and,  by  the  use  of  simple 
functions,  solve  for  the  desired  quantities  in  each  region 
(Ref  13).  For  continuous  structures,  this  involves  divid¬ 
ing  the  system  into  a  finite  number  of  elements  which  are 
interconnected  at  a  discrete  number  of  nodal  points 
their  boundaries.  Displacement  functions  are  chosen  which 
satisfy  the  elemental  boundary  conditions  and  continuous 
forces  or  stresses  are  modeled  as  concentrated  forces  at 
the  nodes.  By  using  the  equations  of  structural  analysis, 
the  forces,  strains,  stresses,  and  nodal  displacements  for 
each  element  are  calculated.  The  global  solution  is  found 
by  appropriately  combining  the  elemental  solutions. 

Constant  Strain  Triangle 

One  element  that  is  used  in  structural  analysis  is 
the  three-noded  triangle,  which  is  also  called  the  constant 
strain  triangle.  A  typical  triangular  element  along  with 
its  nodal  displacements  is  shown  in  Fig.  A-l.  The  nodes 
are  numbered  counterclockwise  as  i,  j,  and  k. 

The  displacements  within  the  element  are  uniquely 
determined  from  the  six  values  of  nodal  displacement.  For 
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the  constant  strain  triangle,  the  function  which  does  this 
is  linear  (Ref  5)  and  is  given  by: 


u  =  +  o^x  +  a^y 


v  =  a.  +  arx  +  a^y 
4  5 


(A-l ) 


The  constants  are  easily  computed  and  are  given  by 


1 

2A 


a . 
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a  . 
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b.  b. 
i  3 


and 
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2A 


c . 
.  l 


a . 
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c  . 
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a  . 
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XJ 


(a-2; 


b.  b. 
i  3 


c . 
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c  . 
3 


Xj 


(A- 3) 


where  A  represents  the  area  of  the  triangle.  The  coeffi¬ 
cients  a^,  b^,  and  c^  and  the  area  2A  are  given  in  terms  of 
the  nodal  coordinates  x^  and  y^  by 

a .  =  x  .y,  -  x.  y  . 

i  irk  k-*3 


b .  =  y .  -  y. 

l  1 2  k 

c .  =  x,  -  x . 

l  k  3 


2A  = 


x 


(A-4 ) 


a.,  b.,  c.  and  a,  ,  b,  ,  c,  are  obtained  by  a  cyclic  -ermu- 
3  3  3  K  K  K 

cation  of  the  indices  in  Eq.  A-4.  One  finally  obtains: 


u  =  IS'1  (ai+bix+ciy)ui  +  (aj+bjX+c^yJUj  +  (ak+bkX+cky)  Uk] 

v  =  ^[(a.+b.x+c.yjv.  +  (a.+bjX+c  .y)  v.  +  (VV+V0  V 

(A- 5) 
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Eq.  A- 5  may  be  written  in  standard  form  as: 


S  *  (v)  -  S  =  IINi.INj,INmUe 


(A- 6) 


ag  represents  a  listing  of  the  nodal  displacements  for  a 
particular  element  and  the  shape  functions,  N. ,  N . ,  and  N,  , 

1  J  K 

are  given  by 


N 


i 


(ai  +  biX  +  Ciy^ 
2A 


(A- 7 ) 


etc . 

The  strain  at  any  point  within  the  element  may  be 
defined  in  terms  of  the  nodal  displacements  derivatives  by: 


Writing  Eq.  A-8  in  terms  of  the  nodal  displacements 
and  coordinates  and  taking  the  appropriate  derivatives 
results  in: 

e  =  B  ae  (A-9 ) 

where  u  is  the  generalized  nodal  displacement  vector  writ¬ 
ten  as 

aj  =  <u.  V;L  u.  Vj  uk  vk)T  (A-10) 
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and 


B  = 


0 


0 

c  . 
3 


0 


c . 
3 

b. 

3 


0 

Ck 


0 


(A-ll) 


Elastic  Analysis 

For  linear  elastic,  isotropic  materials,  the  relation¬ 
ship  between  the  stresses,  a,  strains,  £,  initial  stresses, 
o^,  and  initial  strains,  is  given  by  (Ref  13) 


(A-12) 


where  D  is  the  elasticity  matrix  containing  the  appropriate 
material  properties.  For  plane  stress  D  is  given  by: 


D 


v 

1 

0 


0 

0 

1-v 

2 


(A-13) 


For  plane  strain,  D  is  given  by 


D 


1  v/l-v  0 


E(l-v)  v 

(1+v) (l-2v)  1-v 

0 


1  0 


1-2 
2 (1-v) 


(A-14 ) 


Note  that  for  plane  stress,  a  =  cr  =  a  =  0,  while 

zz  zx  zy 

for  plane  strain,  c  =  e  =  e  =  0  and 
r  zz  zx  zy 


a  =  a  (a  +  a  ) 
x  x  y 
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(A-15) 


t 


As  stated  previously,  before  any  finite  element  analy¬ 
sis  can  be  performed,  the  applied  loads,  which  include  the 


applied  stresses  and  any  body  forces,  must  be  modeled  as 
equivalent  nodal  forces.  For  the  applied  stresses, 


and  body  forces , 


(A-16) 


(A-l'7 ) 


the  equivalent  nodal  forces, 


where 


2i 


(A-18) 


(A-19) 


with  components  IL  and  corresponding  to  the  directions 
of  u  and  v  displacements,  are  given  by: 


q 


e 


a  d(vol) 


b  d(vol) 


(A-20) 


The  generalized  equation  for  the  elastic  response  of 
the  structure  is  derived  from  the  principle  of  virtual  work 
(Ref  13)  and  can  be  written  as: 


K  U  =  P  +  Q  (A-21) 

K  is  called  the  elastic  stiffness  matrix,  U  is  the  general¬ 
ized  displacement  vector,  P  represents  the  externally  applied 
load  vector,  and  Q  is  the  force  vector  caused  by  initial 
stresses  and  initial  strains.  Q  is  used  both  for  time-inde¬ 
pendent  and  time-dependent  problems. 

The  stiffness  matrix  can  be  represented  as: 


K 


D  B  d (vol) 


(A-22) 


v 

The  global  stiffness  matrix  is  obtained  by  evaluating 
the  stiffness  matrix  for  each  element  and  summing  over  the 
whole  structure.  The  equivalent  nodal  forces  due  to  the 
initial  stresses  are  given  by 


Qa 


o 


a ^  d(vol) 


(A-23) 


The  equivalent  nodal  forces  due  to  the  initial  strains 
may  be  expressed  as: 


Qe 


o 


°  ^o 


d (vol) 


(A-24) 


Note  that  the  equivalent  nodal  forces  due  to  the 
applied  stresses,  which  is  the  first  integral  in  Eq.  A-20, 
can  be  written  as: 

Q  _  =  K  ae  (A-25) 

a  ~ 

Once  the  nodal  displacements  have  been  evaluated,  the 
stresses  at  any  point  of  the  element  can  be  found  by  com¬ 
bining  Eqs.  A-9  and  A-12  to  obtain: 
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a  =  DBa  -  D  e  -a 
—  —  —  —  =  —o  o 


(A-26 ) 


where  a  is  the  nodal  displacement  vector.  The  stress 
recovery  matrix  is  defined  as: 


=  D  B 


(A-27) 


Biquadratic  Quadrilateral 

It  has  been  shown  that  for  a  three-no’ded  triangle,  a 
linear  interpolation  function  for  the  displacements  is 
chosen.  Since  the  strains  are  a  function  of  the  first 
derivative  of  the  displacements ,  the  strains  are  constant 
throughout  the  triangle.  This  is  wh.  the  three-noded  ele¬ 
ment  is  called  a  constant  strain  triangle.  The  approxima¬ 
tion  of  constant  strain  in  an  element  is  valid  if  and  only 
if  the  element  is  relatively  small.  For  larger  elements, 
a  higher-order  interpolation  function  must  be  selected. 

Consider  the  distorted  eight-noded  element  shown  in 
Fig.  A- 2 .  While  the  constant  strain  triangle  had  nodes 
only  at  the  vertices,  this  element  has  one  node  in  the 
middle  of  each  side.  The  interpolation  function  chosen  to 
model  the  element  displacements  is  second  order,  and  this 
allows  the  element  to  have  curved  sides.  The  same  interpo¬ 
lation  function  that  models  the  displacements  is  also 
chosen  to  model  the  shape  of  the  element.  Elements  for 
which  this  is  done  belong  to  the  family  of  isoparametric 
elements  (Ref  5) .  The  interpolation  function  may  also  be 
referred  to  as  a  shape  function. 


m 


To  thoroughly  describe  an  isoparametric,  a  nonortho- 
gonal  local  coordinate  system  is  drawn.  This  has  been  done 
to  the  element  in  Fig.  A-2.  The  element  drawn  in  the  local 
system  is  shown  in  Fig.  A-3. 

For  the  equation: 

u  =  N  ae  (A-28) 


the  shape  functions  may  be  given  as  (Ref  12) 

N±  =  i-a+Cq)  (1+nr^)  UCi+nr^-1)  i  =  1,2, 3,4 

=  |(1-C2) (l+nni)  i  =  5,7  (A-29 ) 

=  i(i-n2) (l+Cq)  i  =  6,0 


where  represent  the  coordinates  of  node  i  in  the 

local  system. 

The  shape  of  the  element,  expressed  as  the  coordinate 
points  (x,y) ,  may  be  written  in  terms  of  the  shape  func¬ 
tions  and  the  coordinates  of  the  nodes  as: 


The  displacements  within  the  element,  u  and  v,  may  be 
written  in  terms  of  the  nodal  displacements  as: 
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Figure  A-3 .  8-Noded  Isoparametric  Element 
in  Space 
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In  shortened  notation,  Eqs. 
written  as: 


A- 30  and  A- 31  may  be 


and 


(A- 32 ) 


(A- 3  3 ) 


The  strain-displacement  relation  can  be  expressed  as 


(A-34 ) 


The  derivatives  with  respect  to  x  and  y  must  be  con¬ 
verted  into  derivatives  with  respect  to  £  and  n .  This  is 
done  through  the  use  of  the  Jacobian. 


9-1 


Recall  that  the  Jacobian  is  defined  as  (Ref  12) : 

£ 


J  = 


x,  y , 

n  n 


(A- 3  5 


,  dX  3x 

where  x, „  =  x,  =  5— ,  etc. 

c,  3  c,  n  dn 


The  determinant  of  the  Jacobian  |j|,  is  then  given  by 

'£ 


|j|  = 


x,  y , 

n  n 


(A-36 


=  x,c  y,n  -  x,n  y,c 

The  inverse  of  the  Jacobian,  J* ,  can  be  expressed  as: 


J*  = 


Jii 

J12 

1 

y'n 

•y'£ 

_J21 

J22_ 

|J| 

.  ~x,n 

x'£. 

(A- 3  7 


By  the  chain  rule, 

'<  ), 


5  = 


'£ 


x,n  y'n 


(  ), 


X 


(  ),' 


=  J 


(A- 3  8 


(  )r. 


where  the  empty  parentheses  denote  a  vacant  operator. 
Therefore, 


(A-39 


By  the  use  of  Eq.  A-39,  Eq.  A-34  can  be  written  as: 


10  0 
0  0  0 
Oil 


Jii 

J* 

12 

0 

0 

J21 

J* 

22 

0 

0 

0 

0 

J11 

Jh 

0 

0 

J21 

J22. 

(A-40 ) 


Differentiating  Eq.  A-33  and  applying  it  to  Eq.  A-40 
results  in: 


J* 

J* 

11 

J 12 

0 

0 

j* , 

j* 

_  21 

22 

0  0 

r*  J* 
'21  22 


o  n,c  \v± 

0  N, 


(A-41 ) 


Differentiating  Eq.  A-32  and  applying  the  results  to 
Eq.  A- 3 7  gives  for  the  inverse  of  the  Jacobian: 


J*  = 


N,  V.  -N, _y . 

x  -'ry4-!  -5-^1 


J  ~N,  x .  N,  _x . 

1  1  —  n— i  -  E,-i 


(A-42) 


Since  £  =  B  a&,  the  matrix  is  the  product  of  the  first  two 
matrices  in  Eq.  A-41.  Using  Eq.  A-42  gives: 


N,  y.  -N,ry. 


0  N ,  ^ 

0  %. 


0  0  -N,  x.  -N ,  -X . 

1J|  -  n-i  -  5-i 

--N'n-i  %Zi  -S.£Zi. 

Matrix  multiplication  results  in 


(A-4  3 ) 
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B  = 


%iLiN,rN,ci:iN,n 

0 


-N,  x . N , r+N , rx . N , 

n  i-  5  '5-i-  n 


-N,nx.N,c+N,5xiN,n 


N, 


['n^i-'5_N,5^i-'n 

By  vector  multiplication  identity. 


N,  y. 
—  ry^i 


Therefore, 


N'nYiN,C 


T  „T 
y  .  N, 
■^i  —  n 


T  T 

*i  N'n  -'n 


Defining  the  vector  Q  as: 

2  *  S'n  n'5  *  N'c  % 

(note  that  Qt=:-Q,  that  is,  Q  is  antisymmetric) 
Therefore, 


T 

Yi  2 


N,  y.N,_  -  N,_y.N, 

nJi  5  5Ji  n 


The  B  matrix  may  now  be  expressed  as: 


B  = 


M 


'  T 

4  2 

o 

T  . 
L~— i  C 


0 

T 

-x.  Q 
T 

Yi  a  . 


The  determinant  of  J  may  be  rewritten  as: 


(A-44) 


(A-45) 


(A-46) 


(A-47 ) 


(A-48 ) 


(A-49) 


=  x,cy,n  -  x,ny,c 


T  T  T  T 

xTN , rN ,  y.  -  X.N7  N, -y . 
i  '5  n  i  i  n  5  i 


(A-50) 


97 


Using  Q,  Eq.  A-50  can  be  taken  to  be: 


|J|  =  ~*I  2  ili 

Let 

T 

*x  -  -Zi  Q 
X.  =  Q 

A.y  X  — 


Therefore, 


(A-51) 


(A-52 ) 


(A-53) 


where 


lJl  =  Ky  Z± 


(A-54) 


The  stiffness  matrix  K  can  be  determined  to  be 


K 


D  B  dxdy 


(A-55) 


where  t  is  the  thickness  of  the  element,  assumed  to  be  a 
constant.  By  the  use  of  the  transformation  rule, 

dxdy  =  jjjdCdri  (A-56) 


the  stiffness  matrix  can  be  expressed  as: 


1  1 


K  =  t 


f  f  bt  d  b 
-1  ‘'-l 


J I d^dn 


(A-57) 
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Let 


D  =  D 


D11 

°12 

D13 

°12 

D22 

D23 

(A-58) 

D13 

D23 

°33. 

tdb 

can  be 

written  as : 

B  D  B  = 


Dllxxxx+D33xyxy 

+D13(xxxy+xyxx> 

D12xyxx+D13xxxx 


D12XxXy+D13Xxxx 

+D23xy+D33xyxx 

D22Xyxy+D33XxXx 


(A-59 ) 


+D23xyxy+D33xxxy  +D23  (xyxx+xxxy} 


It  is  the  quantity  B  DB  that  must  be  integrated  in 

order  to  evaluate  the  stiffness  matrix,  K.  While  B  DB  for 

the  three-noded  triangle  contains  only  constant  terms, 

T 

making  the  integration  for  K  trivial,  B  DB  for  the  eight- 
noded  quadrilateral  is  composed  of  variables  in  £  and  n. 

The  integration  for  K  is  no  longer  a  simple  matter  and  a 
numerical  integration  scheme  must  be  used.  The  method  of 
integration  most  widely  used  is  Gaussian  quadrature. 

The  Gaussian  quadrature  integration  rule  in  one  dimen¬ 
sion  (Ref  13)  states  that  the  integral: 

1 


■  / 


f  (C)d£ 


(A-60 ) 


can  be  evaluated  by: 

1 


f  (C)  dC  =  I  H.fU.) 

i-1  1  1 


(A-61 ) 
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denotes  an  appropriate  "weight"  for  each  value  of 


evaluated.  The  location  of  the  sampling  points,  is 


chosen  such  that  the  error  associated  with  the  numerical 
integration  is  minimized.  The  order  of  the  error  is  2n, 
where  n  is  the  number  of  integration  stations  chosen  (order 
of  the  integration) . 

In  two  dimensions,  the  Gaussian  quadrature  integration 
rule  becomes  (Ref  13) : 


I  = 


■I 


n 


f(C,n)d?dn  = 


-l 


n 
Z 

i=l  j=l 


z  HiH.fUj,ni) 


(A-62 ) 


Thus,  if 


f(5,n)  =  bt  D  B  | J | 


(A-63) 


then 


K 


f (^,n)dCdn 


n  n 

t  z  z  H  h  fU  ,n.) 

i=l  j=l  1  3  3  1 

(A-64) 


The  locations  of  £  and  n  are  chosen  symmetrical  to  the  axes 
system. 

The  order  of  integration  needed  to  accurately  deter¬ 
mine  K  is  determined  by  analyzing  the  order  of  detJ  (Ref  5) 

For  a  plane  quadratic  element  of  constant  thickness  detJ 
3  3 

contains  £  and  n  terms.  Therefore,  only  a  2nd  order 


(2  x  2)  Gauss  rule,  where  the  error  is  on  the  order  of  2n, 
or  4 ,  is  all  that  is  required.  However,  when  one  uses  low 


order  integration  rules,  zero  energy  deformation  modes  may 
develop  (Ref  5) .  Consider  an  eight-noded  rectangle  using 
a  second  order  (2  x  2)  Gauss  rule.  Let  the  element  be 
represented  in  £-n  space  and  let  the  nodal  degrees  of  free¬ 
dom  be  assigned  the  values: 


U2  =  V4  =  u6  =  V8  =  0 

(See  Fig.  A-4  for  a  representation  of  the  element  and  the 
displacements . ) 

At  the  Gauss  points,  £=n=±0.5774,  all  of  the  strains  are 
zero.  Since  the  Gauss  points  represent  the  element  when  its 
stiffness  matrix  is  formed,  the  element  offers  no  resistance 
at  all  to  this  particular  deformation  mode.  This  implies 
that  the  resulting  stiffness  matrix  will  be  singular.  If  a 
zero-energy  mode  appears,  it  is  usually  superimposed  on  the 
deformation  modes  formed  from  actual  nodal  displacements. 

This  particular  zero-energy  mode  is  not  compatible 
with  the  same  mode  in  an  adjacent  element.  The  global 
assembly,  therefore,  will  not  have  this  zero-energy  mode 
and  the  global  stiffness  matrix  will  not  be  singular. 

However,  zero-energy  modes  that  are  compatible  may  be 
formed  (Ref  5) .  These  are  detected  by  computing  the  eigen¬ 
values  of  the  global  stiffness  matrix.  For  a  structure 
free  of  rigid  body  motion,  there  will  be  as  many  zero-energy 
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Figure  A -4.  Zero  Energy  Mode  for 

an  8-Noded  Isoparametric 
Element 


modes  as  there  are  zero  eigenvalues  in  the  global  stiffness 
matrix. 

By  selecting  a  third-order  Gaussian  quadrature,  one 
avoids  all  of  the  zero-energy  modes  by  eliminating  the 
extraneous  strains.  The  JALESE  program  makes  use  of  a 
third-order  Gaussian  quadrature  in  performing  the  integra¬ 
tion  necessary  to  evaluate  K.  Figure  A-5  shows  the  eight- 
noded  element  in  £-n  space  with  its  Gaussian  integration 
stations.  Table  A- I,  which  is  Table  8.1  from  Ref  13,  shows 
the  location  of  the  Gaussian  integration  stations  and  their 
associated  weight  functions  for  each  order,  n. 


Table  A 
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I.  Weight  Coefficients  and  Abscissae 
for  Gaussian  Quadrature  (Ref  13) 
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Appendix  B:  Calculation  of  Equivalent  Nodal  Forces 

The  expression  for  the  calculation  of  the  equivalent 
nodal  forces  due  to  applied  loading  and  body  forces  is: 


q 


e 


a  d(vol) 


NT  b  d (vol) 


(B-l) 


where  qe  denotes  the  equivalent  nodal  force  vector,  B  is 
the  matrix  that  relates  strain  to  nodal  displacements,  N  is 
the  interpolation  function  vector,  b  is  the  nodal  force 
vector,  and  a  are  the  externally  applied  stresses. 

Gravity  and  centrifugal  "force"  are  two  types  of  body 
forces.  There  are  no  body  forces  for  the  problem  in  this 
work.  (Gravitational  effects  are  neglected.)  For  an  ele¬ 
ment  at  the  boundary  subject  to  a  distributed  loading  q, 
virtual  work  considerations  result  in  the  equivalent  nodal 
forces,  Feq,  to  be  expressed  as: 

Feq  =  f  q  NT  ds  (B-2) 

*s 

where  ds  represents  the  incremental  path  length. 

For  the  constant  strain  triangle,  the  calculation  of 
the  equivalent  nodal  forces  becomes  straightforward.  Only 
the  integration  of  constants  is  necessary,  and,  for  a  con¬ 
stant  distributed  load  of  magnitude  q,  the  equivalent  nodal 
forces  become: 


Feq 


=  q  £ 


( B—  3 ) 
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where  £  is  the  length  of  surface  that  the  load  covers. 

This  calculation  is  performed  for  each  element  and  the 
loads  at  the  nodes  are  summed.  Figure  B-l  illustrates  an 
example  of  the  calculation  of  the  equivalent  nodal  force 
for  a  constant  strain  triangle. 

For  an  eight-noded  quadrilateral,  integration  must  be 
performed  to  obtain  the  equivalent  nodal  forces.  Consider 
the  eight-noded  quadrilateral  boundary  element  under  a 
constant  distributed  load,  q,  given  in  Fig.  B-2. 

The  interpolation  function,  N.,  for  the  eight-noded 
quadrilateral  is  given  by: 


N. 

l 

=  i(l+CiO  U+T)^)  (C^i+nn^l) 

i  =  1,2, 3, 4 

N. 

l 

=  \a-K2)  (l+nn^ 

i  =  5,7 

N. 

1 

=  (l-n2) 

i  =  6,8 

At  node  point  3,  (£,n)  =  (1,1)  and 

n3  =  i(iH)  (l+n)  (C+n-l)  (b-5) 

Along  the  upper  surface,  n=l.  Therefore, 

N3  =  i(C2+C)  (B-6) 

At  node  point  4,  (£,n)  =  (-1,1)  and 

n4  =  ^(1-S) (l+n) (-C+n-l)  (B-7 ) 

Along  the  upper  surface,  n=l.  Therefore, 

N4  =  (B-8) 
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At  node  point  7,  (5,n)  =  (0,1).  Along  the  upper  sur¬ 

face,  n=l-  Therefore, 


n7  =  jd-C2)  (l+n) 


=  <!-£' 


(B-9J 


Now 


/■ 


Feq  =  Jq  NT  ds  =  SL  I  q  NT  ds  (B-10) 

•'-l 

where  SL  =  length  of  element.  Substitution  results  in: 


Integration  gives: 


i(S2-0  d£  (B-ll) 

(1-52) 


Therefore,  the  equivalent  nodal  force  is  four  times 
greater  at  the  midside  node  than  at  the  corner  nodes. 

As  in  the  case  for  the  constant  strain  triangle,  the 
equivalent  nodal  forces  are  calculated  for  each  element  and 
the  forces  at  the  nodes  are  summed.  An  example  of  the  cal¬ 
culation  of  equivalent  nodal  forces  is  shown  in  Fig.  B-3. 
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Appendix  C:  Users  Guide  for  JALESE 

The  JALESE  program  is  an  elastic-plastic  finite  ele¬ 
ment  code  for  plane  stress  or  plane  strain  problems.  Usage 
of  the  code  is  very  straightforward.  The  data  cards  are 
punched  up  according  to  the  users  guide  by  Ahmad  and 
Papaspyropoulos  (Ref  1) .  There  are  some  fine  points  to  be 
noted,  however,  in  the  loading  and  execution  of  the  prob¬ 
lems  . 

Minimum  usage  of  core  memory  space  is  of  prime  impor¬ 
tance  in  achieving  fast  turnaround.  In  JALESE,  the  largest 
block  of  common  memory  is  used  by  the  stiffness  matrix. 

The  stiffness  matrix  array  is  called  SK  in  JALESE  and  is 
stored  in  the  block  named  NEW2 .  SK  is  a  two  by  two  array 
and  is  sized  by  SK  (2*no.  of  nodes,  1  +  semi-bandwidth) 
(half-bandwidth) .  Therefore,  it  is  important  that  the 
user  size  the  SK  array  for  his  particular  problem. 

To  minimize  the  size  of  the  stiffness  array,  it  is 
important  to  keep  the  bandwidth  small.  This  is  accomp¬ 
lished  by  numbering  the  nodes  appropriately.  For  a  rec¬ 
tangular  array,  the  best  node  numbering  scheme  is  horizon¬ 
tal  (left  to  right) .  If  a  radial  "cobweb"  array  is  used, 
it  is  best  to.  number  the  nodes  along  the  concentric  curves. 
If  a  saw-tooth  pattern  is  desired,  extreme  care  must  be 
used  in  formulating  the  mesh  in  order  to  avoid  a  large 
semi-bandwidth . 

It  has  not  been  found  worthwhile  to  attempt  further 
reductions  in  core  memory  usage  by  reducing  the  size  of 
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the  blocks  that  contain  the  element  connectivity,  nodal 
location,  uniaxial  stress-strain  curve  definition  points, 
etc.  Much  work  is  required  to  insure  that  the  common 
blocks  remain  the  same  length  from  subroutine  to  subrou¬ 
tine.  Eighty  percent  of  the  memory  is  used  to  store  the 
stiffness  matrix.  This  is  where  one  should  direct  one's 
efforts  to  reduce  core  memory  requirements.. 

Besides  minimizing  the  size  of  the  stiffness  array, 
an  efficient  numbering  scheme  for  the  nodes  that  reduces 
the  semi-bandwidth  will  also  minimize  computer  processing 
time.  The  JALESE  program  uses  Gaussian  elimination  to 
solve  for  the  nodal  displacements.  If  the  semi-bandwidth 
is  reduced,  the  computer  will  have  fewer  equations  to 
solve. 

It  has  been  determined  that  the  core  memory  (number  of 
octal  words)  required  to  load  the  JALESE  program  in  its 
published  form  is  270000.  This  is  far  too  large  to  load 
on  the  CYBER  74  (CSB) ;  however,  it  will  load  on  the  CYBER 
175  (CSA) .  Unfortunately,  the  job  will  not  be  run  until 
the  evening  shift  and  if  the  system  is  overloaded  (which 
it  usually  is) ,  the  job  will  not  be  run  until  the  weekend 
shift.  In  the  current  work,  the  JALESE  program  was  modi¬ 
fied  by  cutting  down  the  size  of  the  SK  array  and  the  pro¬ 
gram  was  loaded  with  a  core  memory  size  of  205100. 

It  was  also  determined  that  200  seconds  of  central 
processor  time  and  75  seconds  of  input-output  time  had  to 
be  specified  on  the  job  control  card.  The  exact  amount  of 
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processor  and  input-output  time  actually  used  can  be  deter¬ 
mined  by  examining  the  day file  at  the  end  of  the  print-out. 
One  could,  therefore,  reduce  the  time  specifications  for 
repeat  runs  and  cut  down  turnaround  time.  It  is  advisable, 
however,  to  overspecify  by  25%;  if  the  computer  is  over¬ 
loaded  and  the  job  is  in  the  queues  for  a  long  time,  more 
execution  and  10  time  will  be  used. 

If  a  particular  problem  has  more  than  sixty  elements, 
the  stress  and  strain  tolerances  must  be  increased  from  the 
default  value  of  0.001.  If  the  crack-tip  singularity  ele¬ 
ments  are  used,  one  must  alter  the  tolerance  if  more  than 
fifty  elements  are  used.  If  this  is  not  done,  the  initial 
stress  or  strain  solution  technique  will  not  converge  to  a 
solution.  However,  since  the  program  only  prints  out  four 
significant  figures,  changing  the  tolerance  does  not 
greatly  alter  the  results. 

The  JALESE  program  as  it  is  currently  published  does 
not  have  a  plotting  routine.  The  work  with  VISCO,  which 
does  have  a  plotting  routine,  showed  that  a  plot  of  the 
mesh  was  very  useful  in  debugging  the  input  data.  There¬ 
fore,  a  plotting  subroutine  was  written  for  use  with  a 
modified  JALESE  program.  This  subroutine  is  listed  in 
Fig.  C-l.  To  use  this  subroutine,  it  is  necessary  to  add 
a  line  of  code  that  reads  the  variables  SCAL  and  NSUP. 

This  subroutine  is  compatible  with  the  CALCOMP  565  plotter. 
It  can  be  easily  modified  for  use  with  the  CALCOMP  1036  or 
CALCOMP  1038  plotters. 
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Figure  C-l.  Plotting  Subroutine  for  the 
JALESE  Program 
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As  an  example,  the  following  data  cards  for  the  two 
element  uniaxial  test  case  using  IN-100  are  shown.  The 
user  is  urged  to  examine  the  user's  guide  (Ref  1). 


Card  1:  FORMAT  914,  F20.10 

NPRTY=1 

NELEM=2 

NP0IN=4 

NB0UN=2 

NC0NC=2 

NCP0IN=4 

NQPTS=0 

JPATHS=0 

IREST=0 

GTOL=Q .001 

Card  1’  (used  when  plotting  subroutine  has  been  added  to 
prograr  ) 

FREE  FORMAT 

SCAL=11.  (Scaling  Factor  for  CALCOMP  Plotter) 

NSUP-1  (0  suppresses  printout  of  element  numbers) 

Card  2:  FORMAT  3F20.5,  14 

E=26300 . 0 
P=0 . 3 
YST=130 .0 
NSSPT=4 

Card  3:  FORMAT  4F20.5 

SRS (1) =130 . 0 
SRN(1)=0. 00494 
SRS (2 ) =152 . 0 
SRN(2)=0. 00716 

Card  4:  FORMAT  4F20.5 

SRS (3) =164.0 
SRN(3)=0. 01494 
SRS (4) =164 . 0 
SRN(4)=1. 00000 

Card  5:  FORMAT  314 


NF=1 
NB (1) =1 
NB (2) =1 
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Card  6:  FORMAT  314 


NF=2 
NB ( 1 ) =1 
NB  (2) =0 

Card  7:  FORMAT  14,  2F20.5 
NF=3 

U  (1) =0 . 150 
U  ( 2) =0 . 0 

Card  8:  FORMAT  14,  2F20.5 


NF=4 

U (1) =0 . 150 
U (2) =0 . 0 

Cards  9,10,11,12:  FORMAT  14, 


1=1 

X(1)=0.0 
Y (1) =0 . 0 


1=2 

X (2) =0 . 0 
Y ( 2 ) =1 . 0 


1=3 

X (3) =1. 0 
Y (3) =0 . 0 


1=4 

X  (4 )  =1 . 0 
Y  (4 )  =1 . 0 


Cards  13,14:  FORMAT  914,  F10 

NOD (1,1) =1 
NOD (1,2) =3 
NOD (1, 3) =4 
NOD (1,4) =4 
NOD (1 , 5) =0 
NOD (1 , 6) =0 
NOD (1 , 7 ) =0 
NOD (1 , 8 ) =0 
NEP=1 
THICK=0 . 3 


2F10.5 


NOD (2 , 1) =1 
NOD (2, 2) =4 
NOD (2, 3) =4 
NOD (2, 4) =2 
NOD ( 2 , 5 ) =0 
NOD (2 , 6) =0 
NOD ( 2 , 7 ) =0 
NOD (2, 8) =0 
NEP=1 
THICK=0 . 3 

Card  15:  FORMAT  14,  2F20.10 
N=3 

FX=0 . 600 
FY=0 . 0000 

Card  16:  FORMAT  14,  2F10.10 
N=4 

FX=0 . 600 
FY=0. 000 


Cards  15  and  16  are  repeated  until  the  desired  loading 
has  been  achieved.  The  data  cards  are  terminated  with  a 
blank  card. 

To  run  the  program,  the  following  job  control  cards 
were  used.  Note  that  the  version  of  the  program  used  was 
called  JALESELONGEST.  The  cards  to  activate  the  CALCOMP 
plotter  calls  are  also  shown.  For  more  information  on  the 
CALCOMP  plotter  see  the  CALCOMP  Plotter  User's  Guide. 


HDG,CM207000,T200,I0100.  M790098 ,GANS, 55533 
ATTACH , A , JALESELONGEST . 

ATTACH , CCPLOT , CCPLOT56X, ID=LIBRARY , SN=ASD . 
LIBRARY, CCPLOT. 


FTN , I=A , B=DOG . 


MAP , PART . 

DOG, PL=10000 . 


7/8/9 


Data  Cards 

7/8/9 

6/7/8/9 
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